"""
Code for the third assignment for stochastic simulation,
by Florian Speelman & Jannis Teunissen

Please not that in our code node 0,1,2 are the intersections and node 3 is the outside world,
while in our report 0 is the outside world.
"""

from SimPy.Simulation import *
from SimPy.SimPlot import *
from numpy import matrix
from numpy import linalg
from random import Random, expovariate, uniform

#Store al global variables in a class G(lobal)
class G:	
	rnd = Random(3117)										#Used to obtain random numbers
	roadLength = [500., 1000., 500, 1000., 500., 500.]		#The lengths of the 6 roads, can be different,
	#if the three nodes are i,j,k then the road from i to j has number 2*i + (j > k), etcetera
	roadSpeed = [20000 for i in range(6)]				#The speed cars have on each road, can be different
	carSpace = 20.											#Space per car on the road
	roadR = [roadLength[i] / carSpace for i in range(6)]	#Each road has roadR resources
	trafficLight = [Resource(capacity=1,monitored=True) for i in range(3)]	#Create some trafficlights
	tlRate = [1 for i in range(3)]							#The reciprocal avg. traffic light waiting time
	road = [Resource(roadR[i]) for i in range(6)]			#Create some roads
	roadRate = [roadSpeed[i]/roadLength[i] for i in range(6)]	#The reciprocal avg. road waiting time
	inRate = [0.15,0.1,0.05] 								#Rate at which cars come in at node i from outside
	P = [[0, 0.33, 0.33, 0.33],[0.5, 0, 0.5, .0], [0.33, 0.33, 0, 0.33]]		#Matrix with probabilities,
	#element P(i,j) gives the chance of travelling from node i to node j
	
	#Below are some Monitor() instances we use to keep track of things, names are self-explanatory
	totalCarsMonitor = Monitor(name="Cars in system", ylab="cars")
	timeInSystemMonitor = Monitor(name="Waiting time")
	tlPassedMonitor = Monitor(name="Traffic lights passed")
	nnodeMonitor = [Monitor(name="Queue length") for i in range(3)]

def choose2dA(i,P,g):	#This function comes from an example on the Simpy website
    """  return a random choice from a set j = 0..n-1
       with probs held in list of lists P[j] (n by n)
       using row i
       g = random variable
       call:  next = choose2d(i,P,g)
    """
    U = g.random()
    sumP = 0.0
    for j in range(len(P[i])):  # j = 0..n-1
        sumP +=  P[i][j]
        if U < sumP: break
    return(j)

class Car(Process):		#Cars will be generated to fill our roads and use our traffic lights
	totalCars = 0		#To keep track of the number of cars in the system
	carsProcessed = 0 	#To keep track of the total number of cars that have entered the system
	def __init__(self, startNode):		#Initialize
		Process.__init__(self)			#Call constructor of Process class
		self.startNode = startNode		#The car enters at startNode
		self.node = startNode			#self.node keeps track of the current node
		self.timeInSystem = 0.0			#Total time the car is in the system
		self.tlPassed = 0				#Number of traffic lights passed (or intersections)
		self.startTime = 0.0			#This will store the time at which the car is spawned in the system
	def run(self):						#Time to drive!
		Car.totalCars += 1
		Car.carsProcessed += 1
		self.startTime = now()
		G.totalCarsMonitor.observe(Car.totalCars)		#Since totalCars has changed, we observe it
		while self.node != 3:			#Node 0-2 are the traffic lights, node 3 means you leave the system
			G.nnodeMonitor[self.node].observe(len(G.trafficLight[self.node].waitQ) + len(G.trafficLight[self.node].activeQ)+1)
										#The above line observes the number of cars at the current node
			yield request,self,G.trafficLight[self.node]	#We need to use the traffic light
			t = G.rnd.expovariate(G.tlRate[self.node])		#That takes some time t
			yield hold,self,t
			yield release,self,G.trafficLight[self.node]	#We are done with the TL
			G.nnodeMonitor[self.node].observe(len(G.trafficLight[self.node].waitQ) + len(G.trafficLight[self.node].activeQ))
										#The above line observes the number of cars at the current node
			self.tlPassed += 1			#The car has passed a traffic light
			nextNode = choose2dA(self.node,G.P,G.rnd)		#So let's go to another node
			if nextNode != 3:			#If nexnode == 3 then we leave the system, otherwise we use a road
				nextRoad = 2 * self.node + int( nextNode > (2*nextNode - self.node)%3 )
										#From node i you can take road 2*(i+1) or 2*i + 1
				yield request,self,G.road[nextRoad]			#We need to use the road
				yield hold,self,(1./G.roadRate[nextRoad])
				yield release,self,G.road[nextRoad]			#We are done with the road
			self.node = nextNode							#We have arrived at the next node
		Car.totalCars -= 1									#The car has left the system
		self.timeInSystem = now() - self.startTime			
		G.tlPassedMonitor.observe(self.tlPassed)
		G.timeInSystemMonitor.observe(self.timeInSystem)
		G.totalCarsMonitor.observe(Car.totalCars)			#And after monitoring some more, the car is gone!


class CarGenerator(Process):
	""" generates a sequence of cars """
	def __init__(self, startNode, rate, maxT):
		Process.__init__(self)
		self.startNode = startNode		#The node at which we enter the system
		self.rate = rate				#The rate at which we enter
		self.maxT = maxT				#Until maxT traffic arrives

	def execute(self):
		while now() < self.maxT:
			p = Car(self.startNode)		#We have a new car
			activate(p,p.run())			#Activate the car, so that it's going to drive
			yield hold,self,G.rnd.expovariate(self.rate)	#Wait some time before creating a new car
			
class ResetMonitors(Process):
	""" Resets monitors after some time, to get rid of startup effects """
	def __init__(self):
		Process.__init__(self)
	
	def execute(self):
		G.totalCarsMonitor.reset()
		G.timeInSystemMonitor.reset()
		G.tlPassedMonitor.reset()
		G.nnodeMonitor[0].reset()
		G.nnodeMonitor[1].reset()
		G.nnodeMonitor[2].reset()
		print "Monitors have been reset at time %.2f"%(now())
		yield passivate,self			#The function is done
				

def jackson(tlRate, inRate):
	A = matrix( [	[1., -G.P[1][0], -G.P[2][0]],\
					[-G.P[0][1], 1., -G.P[2][1]],\
					[-G.P[0][2], -G.P[1][2], 1.]])
					#The above matrix satisfies the equation: A * labda[] = inRate[]
	labdaVector = linalg.solve(A, inRate)	#We use a numpy function to solve the linear system				
	
	n = [0., 0., 0.]						#n[i] contains the predicted average number of cars at node i
	maxiter = 10							#Maximum number of iterations to do in the summations
	rho = [labdaVector[i]/tlRate[i] for i in range(3)]	#rho[i] is the utilization at node i
	q = (1.-rho[0])*(1.-rho[1])*(1.-rho[2])	#To speed up the calculation a bit
		
	for i in range(maxiter):				#See our report to understand this formula
		for j in range(maxiter):
			for k in range(maxiter):
				p = q * rho[0]**i * rho[1]**j * rho[2]**k
				n[0] += p * i
				n[1] += p * j
				n[2] += p * k
				
	print "Jackson network prediction of the above:\t%5.2f - %5.2f - %5.2f"%(n[0], n[1], n[2])
	#Now determine the avg. number of cars on each road, and sum
	temp = [(labdaVector[i] * G.P[i][(i+1)%3] / G.roadRate[2*i + int((i+1)%3 > (i+2)%3)]) for i in range(3)]
	meannCars = sum(temp)
	temp = [(labdaVector[i] * G.P[i][(i+2)%3] / G.roadRate[2*i + int((i+1)%3 < (i+2)%3)]) for i in range(3)]
	meannCars += sum(temp)
	meannCars += sum(n)						#Add the number of cars at each intersection to get the total avg. no of cars
	print "Predicted mean number of cars in system:\t%8.2f"%(meannCars)
	print "Predicted mean time in system:\t\t\t%8.2f"%(meannCars/sum(G.inRate))
	

def main():
	initialize()									#Required
	maxTime = 20000									#Run a simulation until maxtime
	resetMonitorTime = 500							#Reset monitors at this time
	
	g = [CarGenerator(i, G.inRate[i], maxTime) for i in range(3)]	#We have three generators, one for each entry node
	for i in range(3):
		activate(g[i], g[i].execute())				#Activate the generators
	
	r = ResetMonitors()								#Reset the monitors after some time to get rid of startup behaviour
	activate(r, r.execute(),at=resetMonitorTime)
	
	simulate(until=maxTime)							#Start the simulation
	
	#Now some output
	print "Total cars processed \t\t\t\t%8d"%(Car.carsProcessed)
	print "Mean number of cars in the system \t\t%8.2f"%(G.totalCarsMonitor.timeAverage(),)
	print "Mean time in the system \t\t\t%8.2f"%(G.timeInSystemMonitor.mean())
	
	print "Average number of cars at an intersection:\t%5.2f - %5.2f - %5.2f"%\
		(G.nnodeMonitor[0].timeAverage(),G.nnodeMonitor[1].timeAverage(),G.nnodeMonitor[2].timeAverage())
	jackson(G.tlRate, G.inRate)

	doPlot = 1										#Show some graphical output if doPlot != 0
	if doPlot:
		plt=SimPlot()
		plt.plotStep(G.totalCarsMonitor,color='blue')
		wtHist = G.timeInSystemMonitor.histogram(nbins = 30, low = 0.0, high = 200.)
		tlHist = G.tlPassedMonitor.histogram(nbins = 10, low = 0.5, high = 10.5)
		plt.plotHistogram(wtHist, title="Time in system", ylab="cars", xlab="time")
		plt.plotHistogram(tlHist, title="Traffic lights passed", ylab="cars", xlab="number of traffic lights")
		plt.mainloop()
	else:
		print "Graphical output surpressed by doPlot = 0"
		
	print "OK"										#Goodbye!
if __name__ == '__main__': main()
	


